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1 Introduction 



When is a gas ionized? 

• Ionization can come from the plasma itself, if hot enough. With X = Ne/Nneutral, 
Saha equation gives, 

T=X = ^(^k B T)^e-^ T , (1) 

where / is the Ionization energy. Comes from statistical physics inside atom + Maxwell 
distribution outside. X ->■ for k B T « /, and X ->■ 1 for k B T > I. 




Figure 1: Degree for Ionization for Hydrogen, with, I = 13.6 eV. 

• Ionization can come from external medium (Ionosphere ? T = say 1000 K). 

• Ionization can come from the proximity of atoms ? Share electrons : metal. 

Classification 

Say temperature T, density N. 

Typical distance between two electrons: iV _1//3 . 

Typical Coulomb energy: e 2 /N~ 1 ^ 3 . 

Typical kinetic energy in classical regime: k B T. 

More kinetic energy than Coulomb: k B T > e 2 /N~ 1 ^ 3 . Big frontier. 

Classical relativistic: k B T > mc 2 . 

Then come quantum effects. When T > Fermi temperature T F , with k B T F = ti 2 (3ir 2 N) 2 / 3 / 2m e . 
Thus, for T < T F , energy increases with density, not temperature. 
k B T F < mc 2 , k B scales like N 2 / 3 . 

k B T F > mc 2 , k B scales like N 1 ^ 3 (White Dwarfs to Neutron Stars?). 
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Important quantities 

Time it takes to neutralize charge in-balance: Plasma frequency 

AirNe 2 



ujI = = 9000 ^/N [cm- 3 l 



p m f 



That's why some waves bounce against the ionosphere. 
Distance over which charge in-balance can exist: Debye length 



Xd = — = \ = 7 A3 10V? 1 [K]/N [cm-3] [cm] 

u p V 47riVe 




Figure 2: Classification of plasmas. 



2 Kinetic theory 



Vlasov and Boltzmann equations 

Say only electrons + fixed positive background. Most basic description level. G(r, p, t)d 3 rd 3 p 
= number of particles in d 3 rd 3 p around (r, p) at time t. How does it evolve? 

Say a particle has momentum p at position r, at time t. Say a force F acts on it. At time 
t + dt, it will have momentum p + Fdt, and position r + p/^m dt (non-quantum treatment). 
Therefore, ALL particles in d 3 rd 3 p around (r, p) at time t MUST be in d 3 rd 3 p around (r + 
p/7m dt, p + Fdt) at time t + dt. That means, 

G ( r + — dt, p + Fdt, t + dt) = G(r, p, t) (1) 

The "hyper" volume element d 3 vd 3 p does not change (Jacobian = 1 here). Just Taylor expand 
the left-hand-side to get the Vlasov Equation, 

dG p dG dG 

G(r,p,t) + - ^dt+ — -Fdt+—dt = G{r,p,t), 

dr 7m dp at 

dG p dG „ dG 
dr 7m op at 

Now, this result is NOT always right. Why? 

We have assumed the force F does not change over dt. But F is an averaged force, in the 
same way the function G is averaged (IGM = 1CT 6 part/cm -3 . If not averaged, d 3 r must be 
gigantic, not infinitesimal). What if something "un-smooth" happened during dt ? 

Close collisions are local, quasi-instantaneous processes, sending some particles OUT of 
d 3 p around p, and some other particles INSIDE d 3 p around p, during dt. (Think about 
billiard ball collisions: local and instantaneous). We'll then have, 

dG p , dG „ , dG , 

G(r,p,t) + —dt+ — -Fdt+— dt 

dr 7m dp dt 

= G(r, p, t) + [Collisions (r, t), q — > p, Vq] — [Collisions(r, t), p— > q, Vq], (3) 
giving the Boltzmann Equation, 

/ [Collisions(r, t), q — > p] — [Collisions(r, t), p — > q]. (4) 

J q 



dG p dG „ dG 
JL_ _| F H 

dr 7m dp dt 



The right-hand-side, referred to as the "collision term", is analitically untractable. Yet, 

2 

that's the one driving the relaxation to a Maxwellian distribution Gm oc e~ v . For practical 
purposes, alternative forms have been worked-out (Fokker-Planck/Landau, Balescu, Krook 
v{G M -G)...). 
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Vlasov or Boltzmann ? 

In a plasma, particles are influenced by, 

• Close collisions, changing p rapidly and appreciably (say 9 > vr/2). Accounted for by 
the collision term in the kinetic equation. 

• "Distant" collisions, which amount to the influence of the overall plasma (p, J — > E, B). 
Accounted for by the Force term in the kinetic equation. 

Define a "close" collision by closest approacb[j] < Rl, such as e 2 /Ri = E K where E K is the 
typical Kinetic energy (fc B T or k B T F ). Frequency for such collisions is roughly v ~ nR^vx, 
with mv 2 K = E K . Time scale for "distant" collisions if ~ u;" 1 . Vlasov's equation, with collision 
term = 0, is valid for v <C u p , i.e., 

n ( f) 2 J** « *» « E K , (5) 

\Ek J V m V m 

which just defines weakly coupled plasmas, where there is more kinetic energy than Coulomb 
potential energy, whether degenerate or not. 



The Vlasov-Maxwell system 

For weakly coupled plasmas, the first equation needed is therefore Vlasov's with F = Lorentz 

dG 



dG p dG 
at 7m or 



E(r,t) + - x B(r,t) 



dp 



0. 



(6) 



System is closed with Maxwell's equations, where charge and current densities are given by, 



p(r,t) = J G(r,p,t)d 3 p, 
J(r,t) = / gG(r,p,t)vrf 3 p. 



(7) 



Eqs. (6j7), together with Maxwell's, form the Vlasov-Maxwell closed system of equations. In 
ID along axis x, we just have for G(x,p,t) and E(x,t), 

dG p dG dG dE A f 

— + — ^ + QE— = 0, —=4nq G(x,p,t)dp, 8 
at 7m dx dp dx J 

with q = — e for electrons. Landau damping comes from these 2, originally with 7 = 1. 



Subscript L stands (again) for Landau. 
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3 From Kinetic to Fluid to MHD Equations 



From Kinetic to Fluid 

Fluid equations can be deduced from the moments of the kinetic equation^ The fluid macro- 
scopic density n, velocity v and pressure tensor P are defined through, 

n(r,t) = /F (r ,u. iM V P(r, t ) = /,*u-v) 8 (u-v)F(r,u, «)-»"», 
n(r,t)v(r,t) = / uF(r, u, t)d 3 u, (1) 



where £g) is "dyadic" product u <g> v = {uivf). If our plasma is cold, which kinetically means 
F(r, u, t) = 5[u — v(r, t)]G(r, t), the density n(r,t) and the velocity v(r,t) are what we 
would expect. Interestingly enough, the pressure tensor vanishes. Microsopic velocity spread 
translates to macroscopic pressure. Consider now the non-relativistic Vlasov kinetic equation, 

dF dF E + v x B/c dF 
at or m av 

The moments of the equation giveQ 

/on d 
[Vlasov] d 3 p — + — .(nv) = 0, 

/q / <9v dv\ / v \ d 

mu [Vlasov] d p =^ mn — + v- — = an E + - X B — — - P. (3) 
\ot or J \ c J or 

For isotropic pressure^] with P = pi, the last term is just the usual gradient dp/dr = Vp. 
The "convective derivative" term (d t + v • V) simply follows a fluid element. 

At this stage, you can close the system introducing a relation between n(r,t) and p(r,t), 
that is, an equation of state. Like for the first moment and the pressure, the Vlasov moment 
J u n ( )d 3 u always yields a macroscopic quantity oc J u n+1 ( )d 3 u from the v • OF /Or term. 
Still regarding the micro/macro duality: a non-zero collision term in the Vlasov equation is 
needed to recover viscosity or friction on the macro level. 

From Fluid to MHD 

We have initially one distribution function Fj(r, u, £) per species. The procedure above shows 
we eventually have one set of fluid equations per species. Assume we just have protons and 

^ee the Appendix of Spitzer's Physics of Fully Ionized Gases for details. Also, Chapter I of William L. 
Kruer, The Physics of Laser Plasma Interactions (Previewed on Google Books). 

2 Not straightforward. See Kruer for details. Note that d/dr is an alternative notation for V. 
3 If the pressure tensor is anisotropic, with P = (pij), 

9 p _ ( dPxX + dp yX ^ 8p zX dp X y + dpyy ^ Op Z y 8p xZ f 8p yZ ( 9j) Z 



dr \ dx dy dz ' dx dy dz ' dx dy 
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electrons of densities n p (r,t) and n e (r,t). If we want to describe fast phenomenon where 
electrons could be decoupled from protons (faster than u~ l , or smaller than A^), we need to 
keep two sets of equations. The so-called Braginskii Equations might be the most elaborate 
version of this option. 

What if we're interested in slow r <C u;" 1 , and large scale 3> Ad, effects? Electrons are 
expected to closely follow protons. The plasma is a electron/proton "soup". Electro neutrality 
on these scales gives n p (r,t) ~ n e (r,t). In the same way we defined the fluid quantities ([TJ) 
and found they obey Eqs. ([3]), we define the MHD variables, 

p{r,t) = m p nJr,t) + m e n e [r,t), V(r,t) = ■ 

m e + m p 

J(r,t) = qn p (r,t)v p (r,t) - qn e (r,t)v e (r,t). (4) 
Combining the fluid equations for electrons and protons yields^} 

dp M 8{pV) 

ar + v 'arj = - xB - v fe + ^) + pg' ( 6 ) 

where pE is neglected with respect to the Lorentz force, as n e ~ n p E ~ 0. Also, a gravity 
term pg is added here. Its fluid counterpart in Eq. ^ would obviously be nmg. The system 
is closed through, 

- = - cV xE, VxB=-J + ^. (7) 

Inserting J = cV x B/47T into Eq. ^ gives the usual magnetic pressure and tension terms. 
The last equation used to close the system is Ohm's law, which simplest version reads 



a 



V 

E + — x B 



where a is the medium conductivity. This equation is just J = crE in the fluid-frame at 
velocity V, transformed in the Lab. frame^J Ideal MHD sets a = oo, so that E = —V x B/c. 
Two concluding remarks: 

• Yes, we sometime consider E = 0, like in Eqs. (|6| and (fright), and sometime E / 
like in Ohm's law or ([7]- left). KulsrucQ explains well how this proves reasonable. 

• We've cheated a little bit. We use the collision/ess Vlasov's equation, and then talk 
about EOS or Ohm's law, which imply collisions. It's just far simpler to forget about 
collisions at the kinetic/micro level, derive the fluid equations, and then get collisions 
back into the game, kind of empirically, at the fluid/macro level. 



4 Eqs. (|3]) formally give a non-linear term n p m p (y p ■ V)v p + n e m e (v e • V)v e 7^ p(V • V)V. An "=" is 
obtained neglecting the electron momentum, and considering V ~ v p . 
5 J.D. Jackson, Classical Electrodynamics, p. 472. 
6 R.M. Kulsrud, Plasma Physics for Astrophysics, p. 44. 
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4 Linear Landau damping - The Maths 

Just a piece of a vast problem: Energy exchange between waves and particles in a plasma. 
Simply put, in terms of the energy transfer direction: 

• Waves — > Particles: Particle acceleration, wave damping. 

• Particles — > Waves: Wave instability. 

The original paper is Ref. pQ. Landau damping is one of the most studied/debated problem 
in plasma physics. Nice Maths and Physical derivatior^j] 

Calculation overview 

Since the calculation is quite subtle and long, it may be useful to get a general overview from 
the very beginning. Here are the steps we will follow: 

1. Derivation of the dispersion equation e(k,u) = from the ID Vlasov-Poisson system. 

2. Landau contour, the continuity requirement and the Laplace transform. 

3. Resolution for small damping and any distribution function. 

4. Maxwellian distribution. 

Dispersion Equation 

Start from ID non-relativistic equations^] for F(x,v =p/m,t) and field E^x^t), 

dF +v dF e EdF (!) 
dt dx m <9t> ' 



dE 

— = 47re 
ox 



n — F(x, v, t)dv 



(2) 



where n = / F dv is the equilibrium density. Assume F = F + F ly with | Fi |<c| F |, F 
being an equilibrium solution. Same for E. The equilibrium electric field E = 0. Linearizing 
Eqs. (l|2), assuming Fi,Ei oc exp(ikx — itut), gives 



= -iuF 1 +ikvF 1 - e— (3) 

m ov 



ikE x 



^ej Fl{x , V , t)dv . (4) 



1 See Kip Thome's Caltech course "Applications of Classical Physics", Chapter 21 mostly for the Maths 
part at http://www.pma.caltech.edu/Courses/phl36/yr2004/- 
2 Easily generalized to 3D. 
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Figure 3: Imaginary part of G = J e u du/u — x, for x G C. The real axis is a discontinuity. 



Extract F\ from the first equation, and plug it into the second, 



e(k, uj) 



0, with, 



e(k, uj) 




k 2 J v — uj/k 



dv, 



(5) 



where up 1 = 47in e 2 /m, f = Fq/tiq and f' Q = dfo/dv. This dispersion relation was first 
obtained by Vlasov in 1925 [2]. It shows uj should be imaginary. Otherwise, we have a 
problem, unless f' Q (uj/k) = 0. The dielectric function e(k,uj) has therefore a real and an 
imaginary part, which for all kind of systems, is related to dissipation. 

We could just consider uj imaginary and take this quadrature as it is, integrating along 
the real axis. But there's a problem. The resulting function of uj is discontinuous, precisely 
when crossing the real axis. As an illustration, Fig. [3] displays the imaginary part of G — 
Je~ u du/ u — x, for x G C. The discontinuity is obvious around Im(x) = 0. One part of 
the plan has to be physically meaningful, and the other not. But which one? We could 
try both options, and check that damping comes only when choosing the upper one. But 
what if we didn't know in the first place that a Maxwellian is stable? We shall see that a 
Laplace analysis of the problem can fully answer the question, and will indeed tell us that the 
"physical" half-plane is the upper one. 

Admitting for now the upper-plane is the physical one, what do we do with the lower 
one? The answer is that we have to "analytically continuate" the function we have on the 
upper-plane, to the lower one. This means finding a function on the lower plane which makes 
a continuous, "analytical" junction, with what we have on the upper one. In this respect, 
a uniqueness theorem from complex analysis helps: if somehow we find an expression in the 
lower plane matching what we have in the upper one, then this is the only one. The Landau 
contour is going to do all of that for us: providing a contour of integration equivalent to an 
integration over the real axis for lm(u) > 0, and an analytical continuation of the later in the 
lower plane Im(cj) < 0. 
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Figure 4: The Landau integration contour. It is not closed. It lies always on the same side of 
the pole. It lies below the pole. 

Landau contour, the continuity requirement and the Laplace trans- 
form 

Let's first give the solution found by Landau, namely the famous "Landau contour". Figure 
[4] shows this integration contour has 3 very distinctive features: 

1. The Landau contour is not closed by the "usual" semi-circle in the lower or upper 
half-plane. 

2. The pole u/k = (u r + i8)/k must always lie on the same side of the Landau contour. 

3. So, which side? The Landau contour goes below the pole. 

These contour prescriptions are called the "Landau prescriptions", and the corresponding 
contour, the "Landau contour". We thus rewrite from now on Eq. (|5j) as 

where J L means integration along the Landau contour. Let's now find out about these 3 
features. 

The contour is not closed 

The contour just goes from v = —00 to +00, and is not closed in the upper or lower half- 
plane, as "usual", because we have no guarantee fb(v) behaves correctly there, so as to cancel 
the integration on the semi-circle at infinite radius. Indeed, considering a Maxwellian with 
/o oc e~ v and setting v = R v e l9v to parameterize the integration on a circle of radius R v , we 
find /o oc e - R v cos2S v which can hardly be considered a vanishing quantity at R v — > 00 for any 
9 V G [0, 7r] (or [0, — 7r], if you close in the lower half-plane). 

Always on the same side 

Assume u = co r + i5, with S > 0. As long as S remains positive in Eq. ([5]), the calculation 
does not pose any conceptual problem as the pole is not on the real axis, and continuity is 
guaranteed. 
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Now, what if 5 approaches 0, and the pole ui/k even gets to cross the real axis? We would 
like e(k, u) to be a continuous function of u. Assume first we leave the integration contour 
unchanged (the real axis for v), and compare the quadrature for to = u r + iS and uj = uj r — id. 
The influence of the pole is mostly felt where the denominator is minimum at v ~ u r /k, so 
let's locally get f out of the integral and compare, 

h= [ , — . rN 7i and J 2 = [ , — , r , , , . (7) 

J v — \u r + io)/K J v — [u r — io)/k 



The difference I\ — I2 is, 



5 Ik 

Jl - /2 = 2z/ FWTWf du 



The continuity of e(k,u) demands the expression above vanishes when 5 — > + . The problem 
is that it does not. Instead, the quadrature tends to n (see function G 2 in Appendix A), so 
that we indeed have a jump of amplitude 2«7r when crossing the real axis^J 

The only way to avoid this is to deform the integration contour in such a way that it 
always lies on the same side of the pole u/k. 

The contour goes below the pole 

To understand why the contour goes below the pole and not above like in Fig. |6j we need 
to follow Landau in rethinking the problem in terms of the time evolution of a perturbation 
applied at t = 0. The Fourier technique is not well suited for that because it entails an 
integration from t = —00 to +00. By design, it does not single out any special moment 
in between. By contrast, the Laplace transform involves times only from zero to +00. As 
shall be checked, the Laplace transform technique gives an unambiguous response about the 
location of the pole with respect to the integration contour. 

Considering a function h(t), its Laplace transform h{u) and the inversion formula 4 ! read 



h{u) = / e tuJt h(t)dt, (10) 
Jo 

h(t) = [ e~ iwt h{u)du, (11) 
Jc L 

where the contour Cl pictured on Fig. [5J passes above all the poles of h(u) at height a > 0, 
and can be closed in the lower half-plane where e~ tU)t behaves conveniently as to cancel the 

3 It can also be said that for I±, the integration path makes a ccwnfer-clockwise half-turn around the pole, 
so that I\ = in. But for I 2 , the half-turn around the pole is clockwise, so that I2 = —in and I± — I2 = 2in. 

4 I here follow Landau's book, [3], p. 139, in defining the Laplace transform this way. That's just the usual 
one, 

9(P)= r 9(t)e- pt dt, (9) 



n 



for p = —iui. It avoids having to rotate everything in the complex plane to relate the calculation to Eq. (|5|). 
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Figure 5: Laplace integration contour. Goes from u = — oo + ia to +oo + ia with a > 0, and 
is closed in the lower half-plane. By design, a > and such that every single poles ui,...,u„ 
of the integrand lie inside the contour. 

integral at infinity there. Note that although the requirement a > is emphasized in the 
book (p. 139), I still have to understand why being above all the poles is not enough. And 
as we shall see very soon, a > is the key to the choice of the right part of the u complex 
plane. 

Let's compute from the Maxwell- Vlasov Eqs. (lj2) the time evolution of the system 
considering, 

F(x,v,t) = riofoW + F^v^e*", (12) 
E{x,t) = E 1 {t)e ikx , (13) 

assuming F%,Ei are first order quantities, and Fi(v,t = 0)e lkx is the perturbation initially 
applied. The linearized Vlasov equation reads, 

9F ^ + ikvF 1 (v,t) - e ^-E 1 {t)f G {v) = 0. (14) 

If we multiply by e tu)t and take the integral from t — to +oo, an integration by part on the 
time derivative term gives, 

f^^^dt = [e^F^t)]"-™ J™ e^F x {v,t)dt 

= -Fx(v,0) -iu}Fx(v,u), (15) 

where lim^oo e tuJt Fi(v,t) = has been assumed. On the one hand, the very existence of the 
Laplace transform of Fi(v,u) = J Fi(v,t)e tujt dt implies it. On the other hand, a important 
conclusion of the paper is that for large times, Fi(v,t) oc e tkvt (see [TJ p. 452, and Plasma 
Talk 5). This point is discussed neither in the book, nor in the original paper. Using Eqs. 
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Figure 6: Forbidden option for the contour. Continuity is preserved, but the contour lies 
above the pole, in contradiction with the Laplace prescription. 



fll5[14| ) then gives. 



en Q 



(ikv - tu)F 1 (v,u) -E l (u) f' (v) = Fi(v, 0), 



m 



(16) 



where Fi(v, 0) now acts like a "source term" at the right-hand-side. A few more manipulations 
exploiting Poisson's equation ([2| give, 



E 1 (co) 



1 4vre r°° F 1 (v,0)dv 
e(k,u) k 2 v-u/k 



(17) 



where e(k,u) is identical to Eq. ([5]). The time dependant electric field given by the inversion 
formula (fTTl) is, 



Ex{t) 



e- i ^ t E 1 (u J )du = 
c L Jc L £\ k M 



-iuit 



4vre f 00 F 1 (v,0)dv 



k 2 



v — ou/k 



duj. 



In contradistinction with Eq. (J5]) where the contour issue is puzzling, the Laplace tech- 
nique used here is clear: The ^-integration in e(k,u) does go along the real axis, and the 



^-integration is performed at fixed Im(w) = cr > 0. It means that in Eq. (18), which com- 
putes a physical quantity, the dielectric function e(k,u) is calculated with u above the real 
v-axis. 

That answers the question we had: the physically meaningful half-plane we were wondering 
about after Eq. ^ is the upper one. The kind of contour pictured on Fig. [6] is thus 
"forbidden". 



Incidentally, what are the poles of the integrand in Eq. (18)? For "normal", smooth initial 



excitations F±(v, 0), the term between brackets won't have poles, so that our poles u}±, . . . , u Tl 



are eventually the zeros of e(k,u). The w-integration of Eq. (18) on the closed contour Cl 
will thus give, with Uj = u r j + iSi, 



Ei{t) = 2iir Res(j) = Aj ex-p(-iujjt) = Aj exp(—iu> r jt)e 
j=i j=\ j=\ 



84 



(19) 



which for large times will be governed by the largest 5j. Therefore, the Laplace transform 
approach cannot spare us the resolution of e(k,u) = 0, as these zeros are the building blocks 
of the temporal response of the system. 
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Resolution for small damping 

We suppose small damping, that is \5\ <C \u) r \, and Taylor expand Eq. ([6]), 



e{k,u r + i5) = e r (k,u r )+iS 



de r 



+ i 



5=0 



€i(k, 0J T ) + io~ 



dei 



8u r 



5=0. 



e r (k, u r ) + i€i(k, u r ) + 5 
de r 



de r dei 



du r doj r 



5=0 



e(k, u r ) + iS 



duj r 



o(6), 



(20) 



5=0 



where the o(5) (negligible with respect to 5), comes from the fact that €i(5 = 0) = (no 
damping, no dissipation, no imaginary dielectric function). 

The first term e(k, co r ) is given by Eq. ^ setting to = co r , or taking the limit of e(k, u r + iS) 
for 5 — > + . The part of the integration along the real axis for v G [— oo, u r /k — e] U [cu r /k + 
e, +oo] gives the so-called "Cauchy Principal Part" , denoted P here. The part corresponding 
to the semi-circle (see Fig. [4] middle) gives the semi- residue for v = to r /k. An alternative way 
of deriving this result, considering the limit 5 — > + , is reported in Appendix A. We thus get, 



e(k,uj r ] 



v 
k 2 



V — LO r /k 



11 dv + mf' {ijj r /k) 



(21) 



This result allows to compute de r /duj r in Eq. (20), which eventually gives, 



e(k, uj) = 1 



k 2 



f'o 



v — uj r /k 

Equating the real part to zero yields, 



dv 



2 r 



k 2 



Trf (u r /k) + 5 



d 



f'o 



v — u) T jk 



dv 



UJ 



k 2 



f'o 



v — uj r /k 



dv = 1, 



(22) 



(23) 



which was the result obtained by Vlasov in the first place. Canceling the imaginary part gives 
directly the damping rate, 

fo(u r /k) 



5 



-7T- 



d p f to 



dv 



(24) 



Eqs. (|23| 24|) formally solve the problem in terms of the distribution function. A first order 

(25) 



evaluation of P~ k 2 /u 2 (see Eq. (27) below), gives 

and then 



UJ r 

5_ 



Up, 



The rate 5 has the sign of f' {uj p /k). That means that if f decreases for v = u p /k, the 
waves is damped because 5 < 0. But if f increases for v = u p /k, we have 5 > and the wave 
can actually grow. 



14 



One could argue we started initially assuming 5 positive, and find it can be negative here. 



It is not a problem for the following reason: Eq. (22) we found assuming 5 > is continuous 
at 5 = 0. It must therefore be identical to the integration on the Landau contour on both 
sides of the real axis. We can therefore confidently solve it regardless of the sign of 5. In 
other words, thanks to the Landau contour, we can compute the result as if 5 was positive, 
and then don't care about the sign. 

Historically, Vlasov first ran into Eq. ([5]). He escaped the problem posed by the pole on 
the real axis by considering only the P of the quadrature. He did so apparently without much 
foundation, which Landau denounced without mercy in pQ. We understand from the analysis 
above that doing so, he missed the imaginary part which would have led to the "Vlasov 
damping" . 



Maxwellian distribution 

Let's finally consider a ID Maxwellian distribution, 

fo(v) 1 



-mv 2 /2k B T 



^k B T/ 



rn 



(26) 



For phase velocities u r /k much larger than the thermal velocity V t h — y/ksT/m, we can 
expand the denominator in powers of kv/u r , since that quantity is small where the numerator 
is relevant. We thus have, 



v — oj r /k 



dv 



, i kv k 2 v 2 k 3 v 3 
f Q 1 + — + — + — 4 



k 2 k B Tk 4 



+ 3- 



+ ... 



to:: 



dv 



(27) 



For small k, namely kV t h/cu r ~ kV t h/uJ p 1, Eq. (23) now gives 

Vk^Tj 



0J„ 



u;(l + 3k 2 X 2 D ), with X D 



LJ.„ 



(2f 



We finally (phew!) use Eq. (24) to extract the damping rate. On the one hand, we compute 
the derivative of the P with respect to u r using Eq. (27), and then simply set u r = co p in the 



result. On the other hand, we set u r = lo p in f to fincj^ 



5 = -u; p ^exp 



2k 2 \ 2 D 



(29) 



Fluid theory just gives the real part of the frequency, namely Eq. (28), so that Landau 
damping is a purely kinetic effect. 



5 Some authors insert the full expression of uj r from Eq. (28), yielding — l/2k 2 \ 2 D — 3/2 in the argument of 
the exponential. 
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G2 

5 




Figure 7: Functions G\ and G\ involved in Eq. (32). For small S/k, G\ is almost 1 everywhere, 
except for v = u r /k, where it is 0. G2 peaks at v = u r /k and tends to elsewhere, while its 
integral is always 71", like a Dirac S function. Parameters are u r /k = 4, 5/k = 0.2 for G\, and 
S/k = 0.2,0.5,1 for G 2 . 



Appendix A 

Let's derive, 



/n -dv = P I ^—dv + nrf' {uj r /k), (30) 



v — uj r jk J v — u r /k 



used for Eq. (21 ), without using the residue theorem. For to = co r + iS with S > 0, integration 
along the Landau contour is equivalent to an integration along the real axis. Let's thus assume 
S > and compute, 

/OO CI 
, .rw, <fo- (31) 

We multiply the numerator and the denominator of the integrand by (v — u r /k) + iS/k, 
which is the complex conjugate of the denominator. We get an expression with a purely real, 
non singular denominator, and clearly separated real and imaginary parts, 

t-v r (v-ur/k) 2 g . r 6 / k w 

" J 00 (v - u r /ky + (6/ky (v - u r /k) dv + 1 (v - u r /kf + (6/ky hdv ~ 

V s ■/ 

V V 

Gi G2 

Regarding the real part, the factor G\ of the integrand is for v = u r /k, and ~ 1 for 
S/k |u — w r /fe|. It tends to the P J f' Q / (v — u r /k) for small S/k (see Fig. [7|. The factor G 2 
of the integrand of the imaginary part departs from only for v ~ u r /k. But its integral is 
always n. For small S/k, the quadrature thus tends to irf^Ur/k), and we are back to (|30)Q 



The limit of i£ with <5 — > + is sometimes written "zO" . The identity 



lim / — — -dx — [ — — -dx — P / — ^—dx + inhia), 

(5->o+ /_„„ x — a — 10 / „ x — a — iv J_ nci x — a 
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This calculation is consistent with the Landau contour integration only for 5 — > + . This 



is because in such case, the real axis along which we perform the integration (32) coincide 



with the Laundau contour. If we were to compute Eq. (32) for 5 — > 0~, we would find the 
opposite imaginary part, just because in this case, the real axis no longer fits the Landau 
contour. The latter, instead, is deformed and keeps passing below the pole, precisely to avoid 
the discontinuity. 
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can be referred to as the "Plemelj Formula" in the literature. For 5 — > , the imaginary part above is 
— iirh(a). 
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5 Landau damping - The Physics, Plasma Echo, and a 
(little) word about the non-linear problem 

While the original paper [I] is purely mathematical, a clearer physical picture is provided 
in Landau's book ([2], §30 p. 126). Suppose we switch on at t = a ID electrostatic wave 
E = E sin(A;x — u;£)x, traveling at v^ = u/k along with a particle with velocity Vq at t = 0. 
For Vq slightly larger than v^, the particle is trapped in the wave potential, where it is going 
to oscillate. Doing so, it ends up with an average velocity close to the wave velocity v^. It 
should thus loose energy, and the energy goes to the wave. Situation is reversed for particles 
initially slightly slower than the wave. They end up gaining energy from the wave. 

If slower particles are more numerous than the faster ones, the wave looses more than 
it gains, which means it is damped. Let's now "Fermi-calculate" this, not following Fermi 
but Jackson [3] and Spitzer [I] (who follows Jackson). Landau uses a slightly different ap- 
proach, still implying a calculation with some small parameter eventually tending to zero. 
I chose Jackson^] precisely because there's no such trick in his strategy. The reasoning is 
non-relativistic. 

To start with, which particles can enter the game? If their velocity is too high relatively 
to the wave, they will flow from one potential crest to another, without much net energy 
exchange. The ones for which energy exchange is possible, are the ones which will be trapped 
by the potential. The wave potential height is, 

Av? = — . (1) 
The maximum particle velocity Av in the wave-frame at v$ must then satisfy, 

-m(Av) 2 = Ap (Avf= 2 ^. (2) 

2 mk 

Thus, only particles with velocity v G [v^ — Av, v^ + Av] in the lab-frame can be caught 
by the wave. This is pictured on Fig. |8j taken from another great work on Landau damping 
by Dawson [5]. For a particle near the center of this interval, we take sinx ~ x for the field, 
and the equation of motion in the wave-frame reads 

d 2 x d 2 x 2 2 qE k 

m ltt? + qE ° kx = ^ -Tp + u osc x = 0, where u osc = ——, (3) 

showing it oscillates in the wave potential with frequency u osc . 

Before it was trapped, the particle energy in the lab-frame was just W\ = mv 2 /2. After 
the trapping, the energy is 

W 2 = \m («J + (t; - v^) 2 ) , (4) 

where \mv 2 ^ is the end translational kinetic energy, and |m(v — v^) 2 can be viewed as an 
internal energy of oscillatior^J If there are nofo{v)dv particles with such velocities, with 



The J.D. Jackson who wrote Classical Electrodynamics. 

2 The energy of a mass oscillating is a potential is split between its kinetic energy and its potential energy 
At the bottom of the potential well, all the energy is kinetic. The term oc (v — v^) 2 in Eq. Q is the kinetic 
energy, in the wave frame, of the particle at the bottom of the well. 
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Figure 8: Division of the plasma between non-resonant and resonant (trapped) particles. Only 
resonant particles contribute to the calculation. After [5]. 

J /o = 1, the energy shift is 

dW = n fo(v)dv(W 2 -W 1 ), (5) 
which we integrate over all particles capable of such exchange, 

AW = / n f (v)dv-m (vj + (v - v^f - v 2 ) . (6) 

Jv^-Av ^ 

Expanding f (v) = f {v^) + (v - v^)f Q {v^) H , the term corresponding to / (^) in Eq. @ 

vanishe^J and we find 

AW^ = -^m W At; 3 /oK)- (7) 



If, like on Fig. |8j particles slower than v$ are more numerous than faster ones, /o(^) < 
and AW 7 > 0, which means particles gain energy at the expense of the wave. The wave is 
damped. We now just have to write that this energy leaves the field E over a time scale a;"*, 

d<yE °dt^ = ~ UoscAW = ^osc^mnov^Av 3 /^^). (8) 
Plugging here the expressions for Av and co osc from Eqs. (2j3 ) we find, 



d(Ei/8n) %V2 ul (El 



dt 3 k 2J{iy 91 V8tt 



28 ('5). (9) 



57T 



3 It does not vanish if you omit the "internal energy" term in Eq. fl4 
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As the field energy oc Eq is damped at 25, the field itself is damped at 5. Setting finally 



uj = u p , we have 



(10) 

identical to Eq. (25) of Plasma Talk 4, up to a numerical pre-factor close to 1 (vr/2 = 1.57 



and 4v2/3 = 1.88). A discussion of the non-Galilean invariance of Eq. (10) is available in 
(p. 180). 

A word on Landau damping and gravitation 

According to Ref. [B], Landau Damping of Gravitational Waves would not be possible. Much 
has been done with respect of Landau Damping of more mundane "gravity waves". The 
stability and vibrations of a gas of stars, by Lynden-Bell, seems to be a quite influential paper 
[7]. The abstract concludes stating "Landau Damping occurs for wave-length smaller than 
the critical one [Jean's]". 

Plasma Echo 

Fascinating consequence of the fact that the density relaxes whereas the distribution function 
does not (many functions have the same integral). Original idea by Gould et al. [To] . 

Suppose we produce an initial electric field perturbation oc e~ lklX in the plasma. The 
Laplace analysis [lj of the distribution function temporal evolution (not the field, nor the 
density) shows it indefinitely oscillates with F = /o + /i(i>) exp(ikivt — ik\x). For large times, 
any velocity integral "phase" -vanishes, 

lim I hiy) e ^t-ihix dv = Q) / u x 

t— >oo J 

which is how we recover zero field and density perturbations. The density perturbation and 
the field associated with fi die out, but f\ doesn't. This is how you reconcile the necessary 
reversibility of the Vlasov-Maxwell system, with the apparent irreversibility of Landau Damp- 
ing. There only seems to be a macroscopic irreversibility, but the evolution in microscopically 
reversible. 

Is it possible to detect this ever oscillating fi(v) at later times ? Yes. Assume we wait for 
a time r, and send another perturbation in the plasma oc e tk2X . The second perturbation is 
going to modulate both / and f\ according to e lk 2~"{t-T)-ik 2 x _ Regarding y 1) we w [\\ recover 
something varying like 



e 



ik\vt— ik\x ^ik2X— ik^v (t— r) ^i(k2 — k\)x+ik2VT~i(k2 — k\)vt (12) 



The key-point here is that contrary to Eq. (11), where k\t 7^ implies the velocity integral 



vanishes at large times, together with the first order density and field, the coefficient of v in 
the exponential above is exactly canceled at time, 
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At this time, the velocity integral will not vanish, and an electric field should reappear in the 
plasma. So you perturb a plasma. You wait until everything apparently calmed down. Then 
you send another perturbation, and at the time prescribed by Eq. (13), an electric field will 
suddenly pop-up "out of nowhere", related to the perturbation you first sent. That is the 
"plasma echo" . 

The idea was experimentally tested soon after the theory came, and the echo was found 
|17j . Mouhot & Villani put it this ways: "A plasma which is apparently back to equilibrium 
after an initial disturbance, will react to a second disturbance in a way that shows that it has 
not forgotten the first one" ([E], p. 40). 

Regarding gravitational systems, Lynden-Bell wrote "A system whose density has achieved 
a steady state will have information about its birth still stored in the peculiar velocities of its 
stars" ([7J, p. 295). 



Nonlinear Landau damping 

We found linear waves are damped. Landau Damping has been experimentally confirmed [8] . 
Here are a few landmarks for large amplitude waves (ID, non-relativistic)|^] 

• Isichenko 1997 [9j: Landau damping valid V amplitude (Theory). 

• Manfredi 1997 pU]: Some large amplitude waves do not decay until t — oo (Numerical). 

• Lancellotti & Doming 1998 [11J: Existence of "critical initial states" for which lim^oo E ^ 
(Theory). 

• Caglioti & Maffei 1998 [12]: Mathematical proof of the existence of some damped 
solutions (Theory). 

• Medvedev et al. 1998 [13]: Damping of waves of finite amplitude and arbitrary shape 
according to e St , with lim^oo 5 = (Theory). 

Mouhot & Villani 2010 [HJ [15] : End of the controversy. Nonlinear Landau damping for 
general interactions, including Coulomb and Newton (therefore also including the case of 
galactic dynamics). 

For any potential V(r) such that |V(k)| = 0(\k\~ 2 ~ e ), with e > 0, and any linearly 
stable distribution function f (x,v), large amplitude perturbations relax in such a way that 
all observables (density, field. . . ), 

*(t) = / f(t,x,v)ip(x,v)dxdv, (14) 

relax exponentially with time. The distribution function itself does not relax to its value at 
t = 0. For small perturbations, f(t,x,v) converges to something that is close to fo(x,v). For 
larger perturbations, the distribution function converges to something that is far from f , or 

4 Thanks to Giovanni Manfredi for the summary! 
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it does not converge at all. The large time behavior of a strongly disturbed solution is still 
an open mystery. 

See [H] for the full report, and a great history of the problem, or [15] for a shorter version. 
Villani was awarded the 2010 Fields Medal for this. 
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6 Beam Plasma Instabilities - Introduction 



Miscellaneous 

From now on, and for a number of Lectures, I'll just go through the Review Paper, "Multidi- 
mensional electron beam-plasma instabilities in the relativistic regime" , Physics of Plasmas, 
17, 120501 (2010). 

Counter-streaming flows, possibly relativistic. Lots of them. Basic system: counter- 
streaming electron beams with nb ,n pQ ,v b0 ,v pQ over a background of fixed protons n^. Main 
hypothesis: 

• Collisionless, Vlasov-Maxwell plasmas (i.e. weakly coupled, see Plasma Talk 2), 

• Homogenous, no boundaries (system size ^> c/u p ), 

• Initially current and charge neutral, = n p ov p o and ribo + n p o = rii, 

• No Bo, to start with. . . 

Motivations: simplest system + Fast Ignition Scenario for Inertial Fusion + Shock Accel- 
eration physics (SNR's, GRB's). See Fig. 2 of Review Paper. 

Particle-In-Cell Simulations: great tool for testing/guiding - See Fig. 3 of Review Paper. 

A multidimensional unstable spectrum 

• 1948: some perturbations with k || to the flow are unstable - Two-stream modes. 

• 1959: some perturbations with k _L to the flow are unstable - Filamentation modes. 
Still 1959: collisionless plasma with T x > T y , unstable for k || y. Weibel. 
Difference between them discussed in Sec. III. F of Review. 

• 1960: some perturbations with k arbitrarily oriented are unstable - Oblique modes. 

Bottom line here: Which one will Nature choose? The fastest. Need to tackle the problem 
globally. 

First: look at flow aligned, then flow-perp and the oblique modes. Second: which one 
grows faster? 
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7 Two-stream Instability 



Two-stream (flow-aligned) modes 

Interesting starting with a cold fluid ID model. Equivalent to Vlasov with fo(v) oc 5(v — vq). 
./Von- relativistic. General case shows it's still relevant for the 3D case. 

Linearize conservation and Euler equations. One set for each electron species, and I omit 
subscripts for clarity. Consider first orders quantities ni p , nit,, E\ oc e lkx ~ wt . Conservation 
and Euler equations read, 

dn d(nv) 

dt + dx ~ ' 1 ' 

dv dv 

m— — \-mv— = qE. (2) 
ot ox 

Once linearized, they respectively give 

ni = n - — — , (3) 



oo — kv ' 
.qEjm 

vi = i r—, 4 



so that, 



Then, from Poisson's equation^ 



qriQ ikEl 

ni = t j — vj. 5 



ikEi = 4nq(rii b + n lp ), (6) 



we get, 

1= U Pb | U PP with ^2 = %Vp i / 7 n 

(u — kv b) 2 (u — kv 0p ) 2 ' p,bp m 

The frequency u — kv is the Doppler shifted frequency. Can't help but thinking it looks 
like an energy conservation equation. Without drifts, v 0t b p = and we just have 

{huf = (hoj pb ) 2 + (hoj pp ) 2 . (8) 

Any ideas? 

Until Eq. @, species are disconnected from each other in the calculation. It is Poisson's 
equation which puts them together, summing the contribution of each species. Assume an 
infinite amount of these, each beamlet going at velocity v, with density n f (v)dv, J f = 1. 
The extension of Eq. Q reads, 

4:irq 2 n fo(v)dv = u%q f f (v)dv ^ 2 = A7iq 2 n 

m(oj — kv) 2 k 2 J (v — oo/k) 21 p0 m 



^^Poisson's equation brings a vectorial equation down to a scalar one. We thus loose information, unless 
k • E = kE. The full 3D analysis shows modes with k± = are precisely like this. 
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identical with the one encountered in Plasma Talk 4, up to an integration by part. 
So if we "Fermi understand" Eq. ([7J), we have everything. 
Introducing the dimensionless variables, 



Eq. ([7]) reads, 



w „ kv ob n ob 
x = — , Z = , a = — , (10) 

OJpp ^pp T^Op 



a 1 , 

+ S . ^o - 11 



{x-Z) 2 (x + aZ) 2 ' 



Diluted beam, a < 1 

There are techniques^] to solve Eq. (11 ) in this regime, always approximately, for all Z. I just 
show here how to find the mode growing the most. 

The beam is just a perturbation to the plasma. The modes of the system should be close 
to the modes of the plasma alone. We thus look for solutions at u ~ u pp , i.e. x = 1 + e. We 
also know that the fastest growing mode should efficiently exchange energy with the beam. 
It should thus have have u/k ~ t>o&. With u ~ u p , that means Z ~ 1. Eq (11) now reads, 

a 1 , 

+ ™ - (12) 



(1 + e-Z) 2 (1 + e + aZ) 2 
As we'll checked, | 1 — Z |<C e and aZ ~aCe, which gives 

a 1 at o a .„„. 

l = ^ + T, l = ^ + l-2e e 3 = -. 13 

e 2 (1 + e) 2 e 2 2 



By setting e = pe %e , we find 



"=(2) • " = -T'°'T' (14 > 



With e ±i2n/3 = -1/2 ± zv^/2, we obtain 3 modes 



= l-^-^" 3 . (15) 



Upp 2 4 /3 2 '- ! 

aV3 



( 16 ) 



= 1 -^73+^ al/3 ' ^tafok. (17) 

As evidenced on Fig. [9j the most unstable mode has Z ~ 1, that is fc/wpp ~ fob- An electron 
from the beam always sees the same electric field. 
How to compute these results in a Fermi-like way? 



1/3 



1/ 



2 See S. A. Bludman, K. M. Watson, and M. N. Rosenbluth, Phys. Fluids 3, 747 (1960). 
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Figure 9: Left: Plot of lm(x) in terms of Z for a = 10 3 , 10 2 , 10 1 and 1. Right: Plot of 
Eq. Q, x 2 = 1 + Z 2 - ^l + AZ 2 . The system is unstable, x 2 < 0, for Z < y/2. 



Symmetric beams, a = 1 



Eq. (11) now reads, 



1 



+ 



1 



{x-Zy ' (x + Z) 2 ' v ' 

which can be solved exactly for all Z, giving 

x 2 = 1 + Z 2 ± y/l + 4Z 2 . (19) 

For Z < v^2, the solution with a minus sign is unstable (see Fig. [9]), with a most unstable 
wave-vector Z m and its frequency x m given by 

Z m = — , x m = + z-. (20) 

For the diluted beam regime, unstable modes are plasma Langmuir waves at u> ~ w p , 
traveling with the beam. Things are not so clear here. The beam is no longer a perturba- 
tion. The waves have Re (a;) = 0, and are the modes of the full counter-streaming system 
"beam+plasma" , each of equal density. 

To wrap-up the most unstable mode characteristics in terms of a G [0, 1] : 



Growth-rate: lm(oo/oj. 



vv) ~ 2 4 / 3 
1 - 



\ (Note that ^| ~ 0.68 > \ 



Frequency: Re(u/u pp ) - ± - ^473 
Most unstable wave-vector: Z = 1 



v/3 
2 ' 



Relativistic effects 

Maxwell's and conservation equations are the same. Euler is now (subscripts omitted), 

d(jv) d(-fv) 
m — — h mv — - — = qhi. 



dt 



dx 



(21) 
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It turns out that when linearizing "7?;" instead of 'V , one finds, 



iv = 7oW + ^i7 H 



As a result, Eq. (11) is replaced by, 

1 



a 



+ 



(x-Zf^ (x + aZ) 2 ^ 



(22) 



(23) 



Intuitively, where does these I/7 3 come from ? If a particle oscillates along its main 
direction of motion, its mass gets a 7 s relativistic boost. Changing m to rwy 3 is Eq. gives 
the result above. 



Diluted beam, a<Cl 

Here, 7 P ~ 1, so that we can recycle the non-relativistic results for diluted beam, formally 
replacing a — > a/75, i- e - n & — * n b/lb- The unstable modes given by Eq. (17) now reads, 



1 a 1 ^ ,V3aV3 



^ = — = 1 " + — • ( 24 ) 



w PP 2 4 /3 7b 2 4 /3 lb 



Symmetric beams, a = 1 



With two symmetric beams, the Lorentz factors are the same 7^ = 7^ = 7. Equation (23) 
now reads, 

1 = (x - Z) 2 7 3 + (x + Z) 2 7 3 ' 

Here again, we just replace x — > X7 3 / 2 and Z — )■ Z7 3 / 2 , and we're formally back to the 
non-relativistic case. Equation pOl) then gives 



Z m -^^, x m -0 + i^-^. (26) 
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8 Filamentation Instability - Part 1 



We still consider the same counter-streaming system, but look now at perturbations with 
k _L to the flow. With respect to the Two-stream instability (k || to the flow), the situation is 
reversed: The physics is simple, but the full maths are involved. Let's start with the physics. 

Physical picture 

Suppose two particle currents of same radius a and density n but opposite velocities u, per- 



fectly overlap (Fig. 10, left). The system is charge and current neutral, in equilibrium. We 



now set them apart by a distance R (Fig. 10, right). The first current generates a B field at 



the level on the second one. The field is such that the Lorentz force F produced repels the 
other current even more. Unstable system. We can write, 

where dM is the mass of the volume element. The force reads, 

7/ 

F = dq-B, (2) 
c 

where dq is the charge of the volume element. With a density n, and particles of charge q 
and rest mass m, the charge dq and the mass dM of the volume dV read respectively, 

dq = qndV, and dM = ■ymndV, (3) 

where 7m is the relativistic mass boost for transverse motion. Equation ([I]) now reads, 

, T ,M n ,,r d2R * M o d2R f»\ 

qndv — B = ^mndv — — r , i.e, q— B = 7m— — -. (4) 
c dt l c dt z 

B is the field created by the current, so that 

21 

B = — , where I = nquna . (5) 
cR 

Replacing the current / by its expression, we find 

— = — , with 6 = u p — =, and f = — , [3 = -, (6) 
dt 2 £ y27 a c 



with = 47inq 2 /m. Although this equation won't give £ oc e St , it does tell the system 
does not relax to its initial state, on a time scale oc which fits exactly the result of the 
linear theory[[] Maybe an exponential grow would be obtained starting from opposite current 
partially overlapping. 

1 \Jp to a factor of order unity, as usual. 
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Figure 10: System unstable to the filamentation instability. 

The Maths 

The Dispersion Equation: Calculation Pattern 

The dispersion equation for the filamentation instability is not easier to derive than the one 
for arbitrarily oriented k's. I will thus go over the general case k\\k± 7^ 0, and then focus on 
fey = 0. For the Two-stream instability, the flow of the calculation wa^j 

Euler + Conservation eqs. (or Vlasov) for each species 

I 

First order densities nj's in terms of Ei 
I 

Merge info for all species through Poisson's equation 

I 

Dispersion Equation 

We could use Poisson's equation for modes with k || to the flow because we know^jthey 
have k || E. In plasma jargon, we say these modes are longitudinal, or electrostatic. For the, 
Poisson equation, which convert a vectorial into a scalar identity, doesn't result in a loss of 
information, precisely because k || E. 

For filamentation modes, we don't know about the respective orientation of k and E. 
The divergence of the electric field would introduce the cosine of the k, E angle, which is 
unknown. Poisson's V • E = 47rp gives k x E x + k y E y + k z E z = Airp, which cannot be used as 
a dispersion equation because it yields only one equation for three components of the field. 
We need a 3D "merging species" equation which does not result in information loss. This is 
Maxwell- Ampere, which merges the currents instead of the charges. 

The general pattern of the calculation is indeed quite similar: 



2 See Plasma Talks 7. 

3 We'll soon find out it is true. 
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Euler + Conservation eqs. (or Vlasov) for each species 

I 

First order currents Ji's in terms of Ei 
I 

Merge info for all species through Maxwell- Ampere equations 

I 

Dispersion Equation 

Let's see this more in details, reasoning again from the fluid equations. Every equilib- 
rium quantities are now slightly perturbed with terms oc exp(?k • r — ioot). The linearized 
conservation equations give for each species: 

ni = n . (7) 

oo — k ■ v 

The linearized non-relativistic (so far) Euler equation give, still for each species: 

-iq/m ( v x B 1 
vi = Jlii H 



oo — k • v 

It is easy to eliminate Bi through Maxwell-Faraday equation, 



Bi = -k x Ei, (9) 

oo 



so that we see how Eqs. (7j8) eventually give n x and vi in terms of E x alone, for each species, 

-iq/m ( v x (k x Ei] 

Vl = : Jiil 



uo — k • v \ oo 

-iqno k / v x (k x Ei) ^ 

n x = : r^- Ei+ . 10 

m {oo — k-v ) z \ oo j 

We may now write Maxwell-Ampere equation, to merge the information from all the 
species into one single equation depending of Ei only, 

— too Aiix 

zkxBi = Ei + — Ji, (11) 

c c 

and eliminate Bi from Maxwell- Faraday Eq. ^ to obtain, 

4i7r 

— kx (kxEi) + Ei + Ji = 0. (12) 

00 z 00 

The first order current is finally expressed through, 

Ji = ^o,feVi ife + n hb v ,b + n 0)P vx, P + nx, p Vo, P ■ (13) 

S v ' S v " 

Beam part Plasma part 
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Although the end result is not really "user friendly" , we can see howEqs. ( ]10ffl2fl3| eventually 
yield a tensorial equation of the form 

T-Ei = 0. (14) 

When starting from the Vlasov equation, linearization gives the first order distribution 
function for each species, 

f /, v n/m ( vxBA df 

/i(k, v, u) = : Ei + ■ — . 15 

<jj — k ■ v \ c / av 

Here again, Maxwell-Faraday Eq. (|9j together with ri\ = J fidv and V\ = J fivdv, allow to 
reach the dispersion equation. 
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9 Filamentation Instability - Part 2 



Dispersion Equation Analysis 

The tensorial equation T • Ei at the end of Plasma Talk 8 has the obvious solution E x = 0. 
Now, the proper modes of our system are precisely the non-trivial solutions T • Ei = 0, with 

That tells us two things: 

• If (3 Ei 7^ / T • Ei = 0) =>■ detT = 0. That's the dispersion equation, yielding u in 
terms of k. 

Assume we pick up one wave vector k. The dispersion equation 

detT(k,w) = 0, (1) 

gives one or more w's, (o>i lk , ■ ■ ■ , ^jv,k) £ Each couple (k, k>j )k ) defines a proper 

mode of the system. Unstable modes have lm(u) < 0. 

The fluid model usually gives a polynomial dispersion equation. Each new ingredient 
to the model (mobile ions, magnetic field,. . . ), adds waves. Polynomial of degree larger 
than 10 are common. 

• The proper modes of the system Ei(k, oS) are in the Kernel of T, which is precisely the 
set of non-zero Ei's fulfilling T • Ei = 0. 

Assume again we picked up one wave vector k. The dispersion equation gives a series 
of frequencies (u^k, • • • , ^jv,k)- We thus have N tensors with vanishing determinants. 
Each of these TV tensors has a Kernel of dimension 1 or 2 (a Kernel of dimension 3 would 
imply T=0). 



T(k,wi, k ) {E M (k,u;i, k )}. =1 



T(k,c^v,k) => {E hi (k,u) N)k )} i=1 



or 2 ' 



or 2 ' 



So, for one couple (k, tUk), the formalism tells how is the Ei field. It lies either along a 
given direction, or in a plane. In particular, the formalism tells us about the k, E angle. 
We don't have to assume waves are longitudinal] (k || E), or transverse (k J_ E). The 
formalism decides for us. 



For a flow || z, and k = (k x , 0, k z ) as pictured on Fig. 11, the final form of the tensor T is 
given by, 

V &xx k z 7] e xz + k z k x 

T= r] 2 e yy -k 2 

rf^xz ~\~ k x k z ff'^zz k 2 

where r\ — u/c and e a/ 3 is given by Eq. (8) of the Review Paper. 



(2) 



1 We could also say we look for the eigen- vectors associated with the eigen-value A = 0. 
2 Also referred to as "electrostatic". 
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X 

k 




»- z 



Figure 11: Axis conventions. 



Two-stream Check 

Let's check our assumption from Plasma Talk 7, that for k || flow, i.e. k x = 0, there are 
longitudinal modes with k || E. Setting k x = in Eq. ^ gives, 



T(k z ,k x = 0) 



V ^xx k 



tfe xz 



ri 2 e xz 

Tj &yy k 

r] 2 e zz 



(3) 



For such wave vectors, the system is perfectly symmetric around the flow axis z. We thus 
have e xx = e yy = e±, ancj^je^ = 0, so that 

r] 2 e± - k 2 



T(k z , k x = 0) 










r] 2 e ± -k 2 




(4) 



The equation T ■ Ei = defines two kinds of waves: 
• Assume (k, u) fulfills, 

e ± = k 2 c 2 /co 2 , (5) 

then, e zz will in general not vanish for the same (k, ui). For these (k, u), the tensor will 
thus have the form, 



0, (6) 
rj 2 e zz ^0 

and waves with with E x 6 (x,y) satisfy T ■ E x = 0. Since k = (0,0, k z ), these are 
transverse modes, k _L E x . In general, they are stable. 



3 Less obvious, but true. 
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If we consider now (k, oj) fulfilling 



£-zz 0, 



(7) 



we find non-zero solutions of T ■ Ei = are waves with Ei G (z), as the tensor now 
takes the form, 

v 2 e ± - k V 

r] 2 e ± - k V 



Since k = (0, 0, k z ), these are longitudinal modes, k || Ei, with dispersion equation, 

which indeed are our two-stream modes. It is thus checked that the modes we investi- 
gated in Plasma Talk 7 do exist. 



The Filament at ion Instability 

About the Dispersion Equation 

Let's now consider k z = in Eq. We find, 







?7 £-xx 

'VV 

^ n ..2.- 12 



r] 2 £ yy - k 2 

rfe zz - k 2 



(9) 



where T • Ei = again defines two kinds of modes: 

• Modes with Ei G (y), therefore transverse since k || x, with dispersion equation, 

k 2 c 2 /u 2 . (10) 



-yy 



The Filamentation modes (at last), with Ei G (x,z) and dispersion equation, 

£xx(&zz k c /oj ) = E xz . 



'ir. 



Of course, we would like to have e xz = 0, which would ease our life and give a simpler, two 
branches dispersion equation, 



o, 

k 2 c 2 /u 2 . 



(12) 



Eq. (12) has been frequently used in the literature to study the Filamentation instability]^] 
It defines purely transverse waves with Ei G (z), that is, || to the flow. The problem is that 
these papers never say they assume e xz = 0. In general, they are wrong. 

I wrote "in general" , because on rare occasions, they study settings for which truly, e xz = 0. 
Which are they? Remember that even if we now focus on k = (k x , 0, 0), this tensor element 



4 See Bret et ai, Phys. Plasmas, 14, 032103 (2007). 
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still depends on the beam and plasma distribution functions. A detailed studjj^] shows e xz 
strictly vanishes only if our counter streaming species are perfectly symmetric. 

So, unless our density ratio is 1, and we have the same temperatures on the beam and the 
plasma, the same Lorentz factors, the same. . . everything, the correct dispersion equation is 



Eq. (11), not (12) 



Cold Analysis - Relativistic effects 

What we've said is so far non-relativistic. Still in the fluid model, the main relativistic effect 
is displayed when linearizing the Euler equation. The relativistic Euler equation reads, 

dp , „. / vxB\ . 

| + (vV)p^ E + — , p = 7 mv. (13) 

Its two linearized versions are, 

im(k • v — w)vi = q ^Ei H — -J , non — relativistic, 

im(k -Vo — u) (70V1 + 7q Vl ^ 2 V ° v j = q + V ° > relativistic. (14) 

Everything is in the anisotropic linearization of 70V around vo- We see above that for a small 
motion along the flow, the relativistic mass increase goes like 7$. But for small motion normal 
to the flow, Vi • v = and the mass increase only goes with 7. This of course, adds a level 
of complexity to the general calculation, as Eq. (10) from Plasma Talk 8 for vi is even more 
involved. 

For the filamentation instability, we have Vi • Vo = 0, and we find we can just formally 
replace m — > 7m. Assuming a cold beam with density rib, Lorentz factor 7^, and cold plasma 
electrons with density n p and Lorentz factor 7p , the tensor elements ara^J 

a 1 

&r.T. J- 



x 2 7 fe x 2 7 P ' 



a 

1 



a aZ 2 1 a 2 Z 2 

ZZ 1 9 7 n 7 • 

x z -f b x^ b x % x 4 7 P 

aZ ( 1 1 



^ 3 7p \ly lb 



(15) 



with again, 



u 7 kvb Hb na\ 

x = , Z = , a = — . (lb) 



ojpp rip 



5 Ibid. 

6 Ibid. 
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Figure 12: Filamentation instability growth-rate for density ratios a — 1, 0.5 and 0.1, from 
higher to lower curves respectively. The beam Lorentz factor is 7b = 10. 



The numerical resolution of Eq. ( jTTj ), when plugging the tensor elements above, yields 
the growth-rate curves pictured on Fig. |l2| As evidenced, the growth-rate just saturates for 
large Z. A trick to recover the large Z growth-rate, consists in extracting the coefficient a n 
of Z n in the polynomial dispersion equation, as a n = is the asymptotic dispersion equation 
for Z — > oo. Doing so, one finds a zero real frequency and 



lim 



8 



Vb a 

C V lb 
c V lb 



, «< 1, 



a 



(17) 



where the agreement with Eq. (6) of Plasma Talk 8 can be checked. Note that for a 



the tensor elements (15) simplify substantially. Equation (12) for unstable modes is valid and 
reads, 

, 2 2ZI Zl 

/3 2 ' 



x 



x 2 lb 



(19) 



which can be solved exactly. 

We can follow the same line of reasoning for the "wrong" transverse dispersion equation 



(12), in order to check its inaccuracy. The exact result is for any a, 



nm 5t =p J*(*lb + l P ) 

Z->00 OJr, 



J pp 



Iblp 




, / «(«7b + l) n „ , 
p\ , a<l, 

lb 



a = 1. 



(20) 
(21) 
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As expected, the result for the symmetric case a — 1 is the same. But for the diluted beam 
regime a Cl, the "transverse" growth-rate 5t differs from the exact one by, 

5 T = V 1 + "76, (22) 

so that the transverse calculation overestimates the growth-rate by a factor which can be 
arbitrarily large. 



37 



10 Oblique Modes 




Figure 13: Axis conventions and setup. 



The electrostatic approximation 

We now come to these fast growing "oblique" modes in the relativistic regime. They are 
found for both k\\ 7^ and k± ^ 0, thus their name. 

Let's remind the general dispersion equation for a set-up like to one pictured on Fig 
From Eqs. (1,2) of Plasma Talk 9, we have the dispersion equation^ 

detT(k,w) = 0, 



13 



with, 





t\ &xz k x k z 



V 2 £xz 



7] Syy ' 





tfe zz 



+ k z k x 


-kl 



(2) 



where r) = u/c and e a p is given by Eq. (8) of the Review Paper. We could just go on with 
this expression, plug some distribution functions for the beam and the plasma, and solve the 
dispersion equation. Doing so, we would realize something important for these oblique modes: 
unless we're really close the ku = 0, these modes have k x E ~ 0. That had already been 
noticed long ago in the first papers on the topiq^J for the cold case. That's been confirmed 
recently for the hot case with various distribution function^ but to my knowledge, it hasn't 
been proved from the formalism. 

It is then fruitful to assume k x E = 0. Although this approximation breaks down for k 
near the normal direction, it has so far been found valid for the fastest growing oblique mode. 
The approximation is called the "electrostatic" or "longitudinal" approximation. 



lr There are some general symmetry requirements on the distribution function. See Review Paper. 

2 Bludman et al, Phys. Fluids 3, 741 & 747 (1960). Fainberg et al, Sov. Phys. JETP 30, 528 (1970). 

3 Bret et al, Phys. Rev. E 70, 046401 (2004) & Phys. Rev. E 81, 036402 (2010). 
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Poisson's equation can still deliver a dispersion equation, but in a slighter intricate way 
because of relativistic effects (a derivation from Ampere's equation, similar to the filamenta- 
tion one, is exposed in the Appendix). We simply go through the calculations of Plasma Talk 
8 & 9, assuming at each steps k x Ei =0, implying Bi = as well. 

The relativistic linearized conservation and Euler equations give for each species, 

k • vi 

ni = n 1 

oo — k • v 

Vl • v \ /_ , v x 



c 2 



m 7oVi+7^v = g E 1 + ^ . (3) 



Solving these two equations gives for the perturbed density. 

.k z E Xz + j 2 k x Et 
(co — k • v 

Inputs from each species are then merged through Poisson's equation 
, ^ 2 k zEu + llKE lx 2 k z E u + lpk x E lx 

K • hi i — UJ b — — r^— , [p) 

p (u-k-vob) 2 pp (w-k-vop)' 
which may be put under the form W • Ei = 0, where the vector W reads 



n l = 1 77. u 777-2 ' ( 4 ) 



Now, k || Ei and W ■ Ei = 0, implies W ■ k = 0, which gives, 

1 _ k l + kill u%hl _ k\ + k 2 xl 2 ul p hl = 

k 2 + k 2 (oo- k z v b ) 2 kl + kl (uj - k z v p f ' 1 ' 

Note that in this longitudinal approximation, there are no oblique effects for 76 = 7 P = 1. 
A generalization of the result to the kinetic level is not as obvious as in the ID theory for the 
two-stream instability, precisely because we are not ID. The kinetic equation reads^J 

= x | W f k -d/o(p)/<9p ^ „ ™ ^ 



k 2 J w-k-v " ^l-v 2 /c 2 ' 

It should not be very difficult to derive intuitively Eq. (|8]) from Eq. ([7]) summing the beamlets 
contributions, as we did for the ID case. In terms of the usual dimensionless variables, 

w kv b n b 

x = — , Z = — , a = — , (9) 



CUpp ^pp 



and using v p = —av b , Eq. ^ reads, 



n _ 1 Z 2 + 1 2 Z 2 ahl Zj+jjZl Ihl 

z} + zi (x-z z y zi + zi (x + z za y { ' 



v v 

beam plasma 



4 The quantity equal to here is called the dispersion "function". See S. Ichimaru, Basic Principles of 
Plasma Physics: A Statistical Approach, Chapter 3. 
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Figure 14: LEFT: Exact growth rate from dispersion equation (18) (transparent) vs. longi- 
tudinal for a = 10~ 2 and 75 = 3. RIGHT: Same for a = 1 and 75 = 1.1. 



Diluted beam 



For a <C 1, % ~ 1 and Eq. (10) reads 

7'2 



= 1 



2 72 

a; 



a/7 



l2" 



'Hi 



This equation is very similar to the one we found for the diluted two-stream case (non- 
relativistic). We formally deal with a diluted beam of equivalent density ratio, 



a 



a 



76 



(12) 



The maximum growth rate will be found for Z 

reads, 



Zl + zi 

1, and the frequency of the unstable mode 



Si 



1 



1/3 



1 + ZJ 

« 1/3 fi + ilzl 



(13) 



2 4/3 



76 



1 + zi 



1/3 



The growth rate (13) displays THE oblique effect: For perp components of the wave vector 



} . s 1/3 

such that 3> 1, we switch from a 7^ to a 7 t scaling. 

The validity of the longitudinal approximation is tested on Fig. 14 where it is compared 
to the exact solution. As expected, it breaks down for small Z z while the exact calculation 
renders the filamentation growth rate as well. 

Fixing Z z = 1, we can compare the exact solution, the longitudinal result and Eq. (13) 



along the perp Z x direction. The result is displayed on Fig. 15 
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Zx 



10 



Figure 15: Comparison of the exact growth rate (top, blue curve), the longitudinal result 



(lower, purple curve) and Eq. (13) [middle, yellow curve] along the perp Z x direction at fixed 
Z z = 1. The density ratio is a = 0.01 and the beam Lorentz factor 75 = 10. 



Symmetric beams 



For a = 1, % = 7 P = 7 and Eq. (10) reads, 

0=1-, 1 



[x 



Z Z ) 2 T^ (x + Z z ) 2 T^ 



7 



Zl + 7 2 Z 2 



1/3 



(14) 



With an equation formally equivalent to the one studied for the two-stream symmetric case 
in Plasma Talk 7. Although the equation above can be exactly solved for x, studying the 
fastest growing Z z for any given Z x is difficult because both are eventually inside T. Once the 
equation is solved, we can however look at the large Z x limit of the growth rate which reads, 



r2 

°Z x ->oo 



l + Z 2 zlb -Jl + AZ^ b 



lb 



(15) 



reaching the extremum, 



2^ 



for 



Vg 
2^ 



(16) 



Exact dispersion equation 



Without the longitudinal approximation, and for arbitrarily oriented k's, we are back to the 
determinant of the tensor ^ for the dispersion equation. It has two branches corresponding 
to the two factors of the determinant, 



-yy 



(r] 2 e xx - k 2 z ){rj 2 e 



k 2 c 2 /u 2 
(j) &xz ~ k x k z ) 



Ei G (y), 
Ei e (x,z). 



(17) 
(18) 
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The second branch therefore holds the two-stream, the oblique and the filamentation insta- 
bilities. As evidenced by the exact plot on Figs. 14 there is a continuous transition from 
two-stream to filamentation modes, probably linked to a common underlying physics. Any 

ideas ? 



Cold hierarchy 

We may finally establish the hierarchy of modes for the cold regime in the (a, 75) phase space. 
The competing modes, with their variation from a <1 to 1, are 

Two -stream ^) = - 1 



2 4 / 3 7. 2 7b 3 



3/2 : 



Oblique, 3(ar) = ^{^j 



1/3 



Filamentation, 5$(x) = — */——>■ —4/ — • (19) 

C V lb C V lb 

The Two-stream case is quickly settled: it is always slower than the oblique unless 75 = 1. In 
the cold regime, the two-stream instability never governs the unstable spectrum^ 

We are thus left comparing oblique and filamentation modes. For the diluted regime, 
the 76 scaling clearly favors the oblique. Situation is more involved near the symmetric 



regime. As evidenced on Fig. 14 RIGHT, the longitudinal approximation gives the good 
order of magnitude for the growth rate, but is not enough to render the "fine structure" of 
the problem. 

What is this "fine structure"? In the diluted beam regime, we clearly have a local ex- 
tremum for oblique vectors corresponding to oblique modes. When approaching the symmetric 



regime, this may no longer be the case. To evidence this, we've plotted on Fig. 16 LEFT 
the growth rate at Z x — > 00 for different parameters. For a = 0.3, we clearly find an oblique 
extremum, which turns to be the dominant mode. But in the symmetric case a — 1, the 
local extremum disappears, giving rise to a monotonous behavior and a system governed by 
filamentation. 

As a consequence, the oblique/filamentation frontier has to be determined numerically 



for large density ratios. The resulting hierarchy plot can be found on Fig. 16 RIGHT. 
The frontier position for a = 1 can be determined exactly from the dispersion equation at 
Z x — > 00. At a = 1, the equation can be solved, and the local oblique extremum vanishes 
when the second derivative of the growth rate at Z z = vanishes. For a ^ 1, this second 
derivative vanishes on the blue line. Note the frontier gets closer to a = 1 for large 7&'s, with 
a convergence numerically found like ^y~°- 395 . 

More detailed in the Review Paper, Section IV-A. 



5 We'll see later that temperature effects change this. This is why the two-stream instability can be observed 
in some real systems. 
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0.2 0.4 0.6 0.8 1.0 1.2 1.4 



Figure 16: LEFT: Growth rate Z x — > oo in terms of the parallel wave vector Z z , for 7;, = 5, 
and a = 0.3 and 1 (lower and upper curves respectively). The local oblique extremum vanishes 
when approaching the symmetric case. RIGHT: Dominant mode in terms of (7^,0). 



Appendix 

We here derive the dispersion equation ([7]) from Maxwell-Ampere equation. The equation we 
used to merge information from each species, namely Eq. (12) from Plasma Talk 8 simplifies 
in the longitudinal approximation, 



x Ei) + Ei 



Am 



00 



-Ji = 0. 



(20) 



From Eqs. (3]4), the first order current is expressed in terms of Ei, leading again to a tensorial 
equation of the form, 

T-E 1 = 0. (21) 
If we have but two counter-streaming species (fream + plasma), the tensor reads^J 

/ / . 1 "pb 1 w lv n n \ 



V 



u + 
-k 



+ 



lb{k z v b -u)) 7p(k z v p -u) 





v b">lb 



+ 



-fb(kzV b -uj) 2 -y P (k z v p -u)) 2 






to (1 








12 

pb 



,2 

°PP 



(22) 



The summing of elements from each species is here obvious again. Note that when k is aligned 
with the main axis, that is k± = or k» = 0, the respective orientation of k and Ei is easily 
determined, because Ei is also found along the very same main axis. Things are here different 
because the orientation of k is arbitrary while the "easy" axis of our tensor are still the main 
ones. 



6 Intcrcstingly, it is not symmetric. I don't understand why. I did check you don't find the correct result if 
you artificially add the missing element [1,3] to make it symmetric. 
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Assume we have Ei fulfilling Eq. (21) and Ei || k. Because T is a linear operator, that 
implies T • k is also the null vector: 

T ■ k = 0. (23) 
The scalar product k • (T • k) must therefore also vanish, 

k • (T • k) = 0. (24) 

The advantage is that the left-hand-side of Eq. ( pij ) is now a scalar, giving us the dispersion 
equation for longitudinal waves with arbitrarily oriented k's. That quantity can be calculated 
from (22), and gives the dispersion equation (J7|. 
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11 Temperature Effects 



I will quickly go through the main temperature, i.e. energy spread effects, on our insta- 
bilities. Let's first start finding out about the limits of the cold regime. 

When are we no longer "cold"? 

The instability process is a matter of wave-particle interaction^} Assume a mode k is ex- 
changing energy with a group of particles. If during one growth period, all particles remain 
in phase with the wave, the interaction is virtually cold. The wave grows as if there was 
no thermal spread at all. Writing that after one growth period, the velocity spread along k 
produces a spatial spread smaller than the wavelength, we find the condition for the validity 
of the cold modeQ 

Av k r 1 < AT 1 , (1) 

where 5 is the growth rate. Note worthily, the condition is not homogenous throughout the 
k space. The spread Av k and the growth rate both depend on k. A given system may be 
virtually cold for the two-stream instability, and hot for the filamentation. 

The same physical picture allows to understand the main effect of temperature. Thermal 
spread reduces the growth rates, precisely because if condition (jl]) is not fulfilled, the wave 
can exchange energy only with a fraction of the particles involved in the cold regime. 

See Section III.C of the Review. 



Two-stream modes 

Assume a ID beam/plasma system with a velocity distribution such as the one pictured on 



Fig. 17 LEFT, and define the temperature parameter, 



For p = 0, we have two counter-streaming symmetric beams. But it is obvious that if p — 1, 
the two distributions make contact, and we end up with a total distribution equivalent to an 
homogenous stable plasma at rest, with velocity spread equal to ±2V. 

Indeed, the dispersion equation is easily computed, and reads in terms of the usual di- 
mensionless parameters, 

1 ~ {x + zy - (pzy ~ {x - zy - { P zy = °' (3) 



A plot of the growth rate is pictured on Fig. [17] RIGHT, evidencing the progressive stabiliza- 
tion of the system for p approaching unity. 

The same pattern holds for more realistic distribution functions. The so-called "Penrose 
Criterion"^] states that distribution functions are unstable if they have more than one local 



1 Though there could be some issues here. See the end of the two-stream section. 

2 Fainberg et al, Sov. Phys. JETP 30, 528 (1970). 

3 Oliver Penrose, brother of Roger Penrose, Phys. Fluids 3, 258 (1960). 
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Im(x ) 




-V V 1 2 3 4 5 

Figure 17: LEFT: Simple toy model for the stabilization of the ID two-stream instability. 
RIGHT: Growth rate in terms of Z = kV/u p for p = 0, 0.9 and 0.999 from higher to lower 
curves respectively. The system is stable for p > 1. 



extremum. Bottom line: for hot enough beam and/or plasma, two-stream can be stabilized, 
relativistic or not0 

Something interesting: It is tempting to relate the former criterion to the formula for 
Landau Damping giving a growth rate oc f^u/k). Nevertheless, we find here unstable waves 
with a distribution function which derivative is almost always zercP] In addition, when the 
system is unstable for p < 1, the real part of the unstable modes is found at u = 0, so that 
f^cu/k) = in our case, while there are no particles at v — 0. This is not an artifact of our 
distribution functions, because u = also with two counter-streaming symmetric Maxwellian 
species. 

I have never seen this kind of issues discussed, except in one single papeij^] There are 
things left to understand. . . 



Filament at ion modes 

Let's extend our toy model consisting in distributions flat up to a certain velocity ("wa- 



terbag"). Consider the 3 distribution functions pictured on Fig. 18 The shaded areas are 
uniformly filled with particles in velocity space. 

• A is a counter-streaming system. Unstable to both two-stream and filamentation insta- 
bilities. 

• In B, we just extend the parallel spread until the distributions come in contact. Ac- 
cording to the previous paragraph, the result is two-stream stable. But the result is 
also anisotropic. Weibel founcQ it is unstable to perturbations with k normal to the 

4 Buschauer, MNRAS, 137 99 (1977). 

5 To be more accurate, the derivative of the distribution function pictured in Fig. 17 LEFT, goes like a 
Dirac' S for v = ±V ±V t . 

6 Phys. Rev. B 43, 14009 (1991), where the problem is pointed out, but not solved. 
7 Weibel, Phys. Rev. Lett, 2, 83 (1959). 
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Figure 18: Distribution functions with different stability properties. 

highest thermal spread: that's filamentation here. B is therefore two-stream stable and 
filamentation unstable. 

• C is built from B, equating the spread in every directions. The result is stable. 

We thus find we can "play" with temperature parameters in order to stabilize some parts of 
the spectrum, or all of it. The bottom line here is that filamentation can be stabilized. 



Note that in Fig. 18 , we carefully tailor the spreads of every species to cancel the instability. 



What if we just play with the beam, that is, the component around v z > on Fig. 18 -A? 
The end result depends on the distribution: 

• For waterbag distributions, filamentation is stabilized beyond a certain amount of beam 
transverse spread. That can be understood as follow^} Assume a current filament of 
radius 1/k and density n^. The current is / ~ qUbk~ 2 Vb- It creates at the surface of 
the filament a field B = Ujck^ 1 . A charge at the surface is therefore pulled in by the 
Lorentz force 

F B = ^ = P)W- (4) 
c V c / 

If there is no temperature, nothing prevents the filament from further pinching, which 
is why the instability extends up to k = oo in the cold regime. But if we're hot, kinetic 
pressure opposes the pinching: A little piece of filament near the surface, with volume 
dV and surface dS, is pulled in by FsnbdV, and pushed out by UbksTdS. Pinching is 
prevented if, 



n^ksTdS > FsribdV 



k > q 



k B T 



x 



n b 



v t b 



(5) 



which is the scaling found from the theory with waterbag distribution^ (we have used 
dV ~ k~ 1 dS ) and v± b is the thermal beam spread.). 



3 Silva et at, Physics of Plasmas, 9, 2458 (2002). 
3 Brct at al, Physical Review E, 72, 016403 (2005). 
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For relativistic Maxwellians, it has been proved that filamentation never vanishes com- 
pletely 10 If Tfe is the beam temperature, the maximum filamentation growth-rate scales 

—77/2 

like T b .In addition, this result still holds for any plasma temperature. 



Oblique modes - general thermal "rules" 

Two-stream modes are unstable up to a finite k m u, and filamentation up to yet another 
finite k m ±. How do we close the unstable domain? Two different answers according to the 
distribution. With waterbag, we close from, and to infinity. See Fig. 10a of the Review. 
With a Maxwellian, these large k oblique modes are stabilized, and we close "normally", as 
pictured for example on Fig. 14 of the Review. 

Oblique modes' temperature sensitivity is intermediate between two-stream and filamen- 
tation. As evidenced in Plasma Talk 10, they tend to be interesting only in the relativistic 
regime, due to the jb scaling which favors them. Two thermal "rules" are useful to grasp the 
influence of the beam spread over the full spectrum: 

• Parallel spread hardly matters. Why? Because in the relativistic regime, it takes a huge 
parallel energy spread to get a reasonable parallel velocity spread. All velocities are 
squeezed against c. See Fig. 13 of the Review. 

• Two-stream does not care about the transverse spread, because particles differing only 
by their transverse velocity will equally stay tuned with a plane wave at k || flow. For 
the same reason, filamentation do care. 

As a result, modes are all the more affected by beam temperature than they are oblique. 

We could rank, from best to worst, the 3 kinds of modes in terms of the way they "resist" 
the various effects: 



Relativistic: Oblique — » Filamentation — > two-stream. 

Density ratio: Filamentation — > Oblique/two-stream. 

Beam temperature: Two-stream — > Oblique — > Filamentation. 

All this ends up with the mode hierarchy pictured on Fig. 20 of the Review Paper. 



The phase velocity diagram - Fig. 17 of the Review 

A great tool to understand the physics. Plot, on the very same graph, the distribution 
functions and the phase velocities of the unstable modes. On can straightforwardly check 
which species are in resonance with which kind of modes. 



Gremillet, Unpublished, Bret et at, Physical Review E, 81, 036402 (2010), 
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12 Non-Linear Regime 





Figure 19: Left: The "Mandelbrot set", defined by Eqs. (Tp). Right: Typical field growth 
extracted from a PIC simulation. Growth roughly stops when the linear exponential growth 
stops. The line shows the expected growth from theory. 

Non-linear problems may be complicated with just 1 or 2 degrees of freedom. Think about 
this 2D example: take the very simple sequence^ 



z n+l 



Z n + C i 

c, c G C, 



and just define the set 



M 



|c6C / lim | z n 1 7^ oo \ . 



(2) 



You have the famous, and incredibly complicated, fractal and everything, "Mandelbrot set" 
pictured on Fig. 19 People were stunned when they realized something as trivial as Eq. (jlj 
could generate such amount of complexity]^} 

A plasma has oo number of degrees of freedom. Yet, to my knowledge, something as 
beautiful as Fig. 



19 is still lacking in plasma physics. Maybe because it's too complicated. 



At any rate, there's no hope of analytically finding out about the long term evolution of 
our beam plasma systems in the general case. Remember it took a Fields Medal to prove 
non-linear Landau damping. Even for the cold case, things are not easy. 

I'll go through some results on the saturation of the various instabilities, always assuming 
the fastest growing mode is the only one excited, and that everything is cold at t = 0. The 



^on-linear in the sense that if two sequences fulfill z n+ i = z 2 + c, their linear combinations do not. 
2 A Math-guy friend of mine once told me people would laugh at Mandelbrot, as "the guy who works on 
polynomial of degree 2" ! 
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former is quite reasonable, as the most unstable mode grows exponentially faster than the 
rest. Relativistic effects help as they drive a sharper unstable spectrum, where growth rates 
vary rapidly from one mode to another. The later is a limitation. 

The idea is to find out when the linear theory should break down, and to claim that 
growth stops at that point. Granted, exponential growth should stop there. But other kind 
of growth could keep on. Yet the observed field growth in PIC simulations is always like to 



one pictured on Fig. 19 Why? 



See Section V of the Review Paper, and references therein. 



Two-stream and oblique instabilities 

Two-stream, non-relativistic 

We assume a diluted beam. There, two-stream is resonant with the beam. Say the wave 
Ee lkx ~ luJt is growing, traveling along with the beam electrons. Electrons will start oscillating 
in the field at the "bouncing" frequenc}|^J 

u b = . 3 

m 

The linear assumption that all electrons have v = v b during one growth period S^ 1 breaks 
down when, 

w> = 6 =► B.= *??*, (4) 
u p q 

giving the value of the field at saturation (I've set here k ~ u p /vb). A great by-product of E s 
is the beam energy loss AWb. Since the field energy can only come from the beam energy, we 
can write, 



AW b E 2 J8it 



\n b mvl \n b mvl 



(5) 



n b \u p ) \n p 



1/3 



where I just replaced the growth-rate by its cold value (%/%,) 1//3tl V We could even write 
the energy lost is shared between the plasma and the fielcQ In such case, AWb/Wb is half 
the result 

As the field grows and traps the beam electrons, they start oscillating in the wave potential. 
PIC people love to plot density graphs in the (x,v x ) phase-space (say x is our dimension) 



such as the one pictured on Fig. 20 At t = 0, the beam and the plasma are just two lines. In 
the instability phase, the plasma does not move a lot, but the beam particles start oscillating 
in the wave, creating "holes" in the phase-space. 



3 See Plasma Talk 5 on Landau damping. 

4 Lorenzo Sironi told me you see this in the PICs. But why? 
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Figure 20: Phase-space hole. 



Two-stream, Relativistic 



A naive reasoning gives the right answer. As the oscillatory trapping motion is along the flow, 
we replace m — > 7^m in the bouncing frequency This gives the field at saturation^ 

E s = 6 -^. (6) 
u p q 

The relative energy loss is then computed as in Eq. ([5]), replacing the growth-rate by its 
relativistic value. It reads, 

^V^(fY /3 . (7) 

Of course, this quantity has to remain small since linear regime means unperturbed trajecto- 
ries, that is, AWb <C Wi,. The energy loss eventually relies on the parameter S witl 



If you do the "clean" calculation going to the wave frame, like in the footnote, you need the 
beam dynamic in that frame to be non-relativistic in order to compute easily a bouncing 



Let's denote quantities in the wave frame with a prime. If the beam dynamic in the wave frame is 
non-relativistic, we have 



12 

0J h = 



qEk' 



m 

where the field E does not need prime as it is parallel to the motion. The bouncing frequency in the lab frame 
is simply = uj^/jf,- In addition, fc = jf,k' since in the wave frame, u/ = 0. We thus find 



qEk 



2 

b = ~3 
76 



rn 



retrieving the m — > j^m "rule". 
6 See Review Paper. 
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frequency. Doing so requires S <C 1, as explained in the Review. For arbitrary S"s, one has 

AW b S 



yielding a maximum energy loss in the linear phase 

AW b 6 /3 „ _ „ .2 



(9) 



.,.-0.18 for S = -. (10) 
VV 6 maa; 25 V 5 3 

Oblique instability 

Poorly known. What is known is that Eq. d9|) roughly works until S ~ 0.45. For larger S"s, 



AWb/Wb seems quite insensitive to S, and remains around 0.18 given by Eq. (10). 



Filament at ion instability 

There are 3 different ways to evaluate the field at saturation! 

1. Cyclotron frequency = growth rate. Filamentation instability grows magnetic 
field. Such field B affects particles on a time scale given by the cyclotron frequency, 

Q Bl ft i ^ 

u c = . (11) 

Writing again that the linear regime keeps on until u c = 5, we get, 

B« = ™^5. (12) 

If one tries to compute the relative energy loss ^/n^mc 2 , the cold result gives a 
factor of order /3 2 , without any more scaling in jb, and almost none in a = nb/n p (see 
Table [D. The conclusion is that estimating the energy loss requires a finer calculation 
than this one, and that the result is quite stable in terms of these variables. 

2. Bouncing frequency = growth rate. Historically, the field at saturation has rather 
been evaluated this way. With v ob = (0,0,i>&) and k = (k x ,0,0), the growing magnetic 
field reads B 1 = (0, B\ sin k x x, 0). At first order, Newton's law projected on the x axis 
gives, 

dj^X Vb 

lbm ~dt? = V B i~ sink x x - ( 13 ) 

Particles at x ~ oscillate at, 

2 = qB^^ 

7&mc 

Note that we assumed there's only a B field here. We know that unless the system is 
strictly symmetric, it's wrong. Here again, we can claim exponential growth keeps on 
while the motion is almost unperturbed, that is until 

qB x v h k x 2 %mc 5 2 

= 6 =>• B s2 = — . (15) 

Ibmc q v b k x 
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Criteria B s B s , mcujp/q units AW&/W& 

Cyclotron frequency = growth-rate ™^b~ j3[a(a + 1)7&] 1/2 /3 2 (a + l)/2 

Bouncing frequency = growth-rate "^y^^^r ,5a (a + 1) /3 2 a(a + l) 2 /27b 

Larmor radius = characteristic k~ x ^^^Vokc f3^ h (3 2 %/2a 

Table 1: Summary of the 3 ways to evaluate the field at saturation B s for the filamentation 
instability, considering k± = u p /c. Results for the cold case. Taking k± oc 7 b , the 3 criteria 
give the same scaling for a = 1 pQ. 



Comparing Eqs. ( 15fl2 ) gives, 

B s2 = B sl - 5 -, (16) 
v b k x 

so that both estimates give the same result only with k x = 5/vf,. For example, with the 
symmetric cold case where 5 = u p f3^2/ r ] h , that implies 

** = (17) 

For non-relativistic setting, the k x is the typical expected one. For 75, very recent cold 
PIC's pQ found indeed that it is the fastest growing k±. Why exactly, as the cold growth 
sate juste saturates at large fcj_'s? 

Larmor radius = characteristic l/k±. Equating the Larmor radius of an electron 
in a field B to the characteristic k c of the insta gives, 

o ^mc v k c 
B s3 = v k c = B s i 



8 v ' 

The 3 results are summarized in Table [Tl considering the cold symmetric case, and taking 
k = 0J p /c for the typical k±. Taking k± oc 7 6 , the 3 criteria give the same scaling for a = 1 



Fate of the filaments 

Opposite filaments repel, but like filaments attract. In our 3D world, filaments turn around 
each other, and like filaments merge. The merging process has been modeled, and successfully 
simulated with PICs3 



More realistic settings, successive instabilities 

We've been so far interested in the short term evolution of the system, that is, the end of the 
linear phase. What's next? I'll just comment Fig. 40 of the Review Paper. The initial setup 
was: 

7 Medvedev et ai, The Astrophysical Journal 618, L75 (2005). 
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Beam: Maxwellian, rib — n p /lO, 7& = 3, = 50 keV. 
Plasma: Maxwellian, T p = 5 keV. 

• < Upt < 80: System initially governed by oblique modes. E field grows at 0.07o;p. 
Heating "kills" de oblique. 

• 80 < u p t < 160: System switches to a two-stream regime. E field grows at 0.016a;p. 
Heating "kills" two-stream. 

• 200 < Upt < 600: Remaining drift feeds filamentation. B field grows at 0.005w p . 

By the end of the simulation u p t ~ 600, the beam had lost about 30% of its energy, 
entirely transferred to plasma electrons. Open questions: 

• Is filamentation the necessary end state of every initial setup? 

• Does the drift eventually ends (in the frame of the fixed ion background)? That is, is 
the drift energy eventually converted at 100% into heat? Sounds reasonable. Is that 
sure? 

Indeed, the interesting question might be how long does it take? 

References 

[1] A. Bret et ai, Submitted to Physics of Plasmas, (2012). 
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13 Ohm's law and the Biermann battery 

The original 1950 paper is Ref. [I], "Uber den Ursprung der Magnetfelder auf Sternen und 
im interstellaren Raum", published in Zeitschrift Naturforschung Teil A. Cited more than 200 
times, and probably read by no one but the happy few who 1) read German and 2) could 
access it. 

■ The MHD equations are formed from the fluid equations for electrons and ions (dropping 
subindices), 

dn d . s 

/ dv <9v\ / v \ 

mn I — — h v • — = gnlEH — xB - Vp + nmg. (1) 



\dt dv ) 
The MHD variables are defined from n ei (r, t), v ei (r, t) a^J 

p(r, t) = miUi + m e n, 
J{r,t) = qniVi-qn e v 



r e, 



V(r,t) = 1 (n e m e v e + m^Vj). 
Merging the fluid Euler equations for both species gives, 

P {~dt +Y ' ~dr~) =- xB - V fe+^)+Pg- ( 2 ) 

Ohm's law is the equation giving the current J. Where does it come from? What we did 
in Plasma Talk 3 was to follow the basic, "business as usual" procedure: sit in the frame of 
the fluid locally at V. There, the electric field is E'. Ohm's law gives the current in the lab 
frame from E', as J' = crE'. Because E' = E + V x B/c, we get the famous (non-relativistic) 

J = a (e + ^- xB^j . (3) 

What about J' = crE'? It arises from the microscopic picture that under the action of 
an electric field, particles, mostly electrons, are accelerated in the direction of the field, while 
collisions with the ions act like a friction force^J Writing something like m<9jv = gE' — w 
where v is some collision frequency, and setting dt = 0, indeed yields v = ^E' and then 
J' = qn^~E'. Ideal MHD assumes v — 0, giving 

E' = E + — x B = 0. (4) 

c 

1 In Plasma Talk 3, the MHD velocity V was defined through (m e + »rii)V = m e v e + mjVj. Most books 
[21 [3] present definition ^ above. It is more rigorous, as it gives a pV MHD term exactly equal to the total 
momentum. At any rate, the difference between the two quantities is oc (n e — rij)(v e — Vj). This is a second 
order quantity in the MHD regime, where electrons are expected to closely follow the ions, so that n e ~ m 
and v e ~ Vj. 

2 Solid state physics call this the "Drude model", from Paul Drude, who came up with this idea in 1900. 
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Figure 21: How a finite size shock can generate non-parallel density and temperature (hence 
pressure) gradients. From |3], p. 406. 

■ But we inadvertently assumed many things. For example, d t = assumes E' varies 
slowly enough with time. If variations are too fast, the stationary regime does not have 
enough time to set in, and some u dependency appears. Also, we assumed particles are 
accelerated only by E' between two collisions. What if B is strong enough to curve the 
trajectories in between? In this case, the resistivity in the direction normal to the field is 
higher (xl.9) than along the field ([2] p. 28, or [3] p. 43). 

Since conductivity comes from the electrons^ let's write their full Euler equation, 

meUe \ W + Ve ' Ijr) = qeUe ( E + ~c X B ) ~~ VPe ~ U ( Ve ~ ^ 

Neglecting the left-hand-side for the small electrons inertia, and setting v = for ideal MHD 
yields, 

E^xB = ^. (6) 

c q e n e 

According to Kulsrud [3] p. 405, the pressure term above is negligible when there is a B field. 
WHY? But for small B's, or even B = 0, you need to keep it. Inserting the electric field 
above in dtQ = — cV x E and setting v e ~ V gives^J 

_ = Vx(VxB)+c^ r ^. (7) 

The second term is our Biermann Battery. There's no B in there, so that it can make it from 
nothing. Still, the microscopic derivation shows we need ionization, just to be able to create 
electronic currents moved by the electronic pressure. 

■ But, an equation of state usually gives p e (n e ). The gradients for p e and n e are thus 



3 The a of J' = ctE' has the mass on the denominator. 

4 We need a little drag between ions and electrons to write v e ~ V. 
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likely to be parallel, so that the cross product vanishes. How do we break this? Several 
possibilities, 

1. Try a system in rotation around z, with Q(z). The pressure gradient reads VP oc 
fl 2 (z)p(r, 9). It has a z component, while Vp(r, #) has not. Those who've read the 
paper say Biermann considered this option. 



2. Try the scheme of Figure [21] Suppose a finite size shock travels through a cold, diluted 
upstream. The downstream just behind the shock is dense and hot. Far downstream 
though, the plasma expands but remains warm, so that it is now warm and diluted. 
The figure shows well how this can generate non-parallel density and temperature (hence 
pressure) gradients. A non-spherical shock also does the job for similar reasons. 

Everything is easily adapted to a system partially ionized. If n n is the density of the 
neutral, define x — n e/ ( n i + n n)- Eq. ^ can be adapted dividing the Battery term by 1 + x 
(0, p. 406). See [7j. 

Kulsrud proposed this mechanism to produce Cosmic Fields from scratch in 1997 [5]. Yet, 
he co-authored in 1992 another paper [6] advocating spontaneous plasmas fluctuations. 

Also Khanna [8] , showed that a rotating BH in a plasma will always generate toroidal and 
poloidal magnetic fields. 

Ref. [9] studied the Biermann Battery effects in Cosmological MHD Simulations of Pop- 
ulation III Star Formation. In its own terms, "We find that the Population III stellar cores 
formed including this effect are both qualitatively and quantitatively similar to those from 
hydrodynamics-only (non-MHD) cosmological simulations" . No dynamical effects. 

The Biermann battery mechanism was successfully tested in the lab in 2012 |10j . 
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